{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Plots.PyPlotBackend()"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "using DifferentialEquations, Plots\n",
    "using Revise\n",
    "using DiffEqOperators\n",
    "pyplot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "telegrapher (generic function with 3 methods)"
      ]
     },
     "execution_count": 62,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Obtained from https://en.wikipedia.org/wiki/Telegrapher%27s_equations#Lossy_transmission_line\n",
    "function telegrapher(L, C, R=0, G=0 ;n=100, vs=t->1.0, is=t->0.0)\n",
    "    if L == 0 || C == 0\n",
    "        error(\"Both the inductance and capacitance need to be non-zero.\")\n",
    "    end\n",
    "    return function(t,u,du)\n",
    "        B2(t,view(u,n+1:2n),view(du,1:n))\n",
    "        scale!(view(du,1:n), -1/L)\n",
    "\n",
    "        B1(t,view(u,1:n),view(du,n+1:2n))\n",
    "        scale!(view(du,n+1:2n), -1/C)\n",
    "        \n",
    "        @. u[1:n] += G*u[1:n]\n",
    "        @. u[n+1:2n] += R*u[n+1:2n]\n",
    "        \n",
    "        # Update the boundary conditions after every iteration according to dependency on time\n",
    "        update_coefficients!(B1, BC=(vs(t), vs(t)))\n",
    "        update_coefficients!(B2, BC=(is(t), is(t)))\n",
    "        \n",
    "    end\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "tline_solve (generic function with 3 methods)"
      ]
     },
     "execution_count": 54,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function tline_solve(u0,vs,is,L,C,R=0,G=0;t=(0.,2.),n=100)\n",
    "    dx = 1/(n-1)\n",
    "    B1 = LinearOperator{Float64}(1,2,dx,n,:Dirichlet,:Dirichlet;BC=(vs(0),vs(0)));\n",
    "    B2 = LinearOperator{Float64}(1,2,dx,n,:Dirichlet,:Dirichlet;BC=(is(0),is(0)));\n",
    "\n",
    "    prob = DifferentialEquations.ODEProblem(telegrapher(L,C,R,G,n=n),u0,t)\n",
    "    solver = CVODE_BDF()\n",
    "    sol = DifferentialEquations.solve(prob, solver, reltol=1e-8, abstol=1e-8,saveat=0.01)\n",
    "    return sol\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [],
   "source": [
    "n = 500;\n",
    "u0 = vcat(1,zeros(2n-1));\n",
    "soln = tline_solve(u0, t->1.0, t->0.0,2,2;n=n);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Base.RefValue{Tuple{Tuple{Float64,Float64,Any},Tuple{Float64,Float64,Any}}}(((1.0, 0.0, 1.0), (1.0, 0.0, 1.0)))"
      ]
     },
     "execution_count": 45,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "B1.boundary_condition"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<img src=\"\" />"
      ]
     },
     "execution_count": 57,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "plot(soln(2.0)[1:n]) # Voltage distribution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<img src=\"\" />"
      ]
     },
     "execution_count": 58,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "plot(soln(1.0)[n+1:2n]) # Current Distribution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<img src=\"\" />"
      ]
     },
     "execution_count": 59,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "plot(soln(0.5)[1:n]) # Voltage distribution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<img src=\"\" />"
      ]
     },
     "execution_count": 60,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "plot(soln(2.0)[1:n]) # Voltage distribution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<img src=\"\" />"
      ]
     },
     "execution_count": 61,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "plot(soln(0.5)[n+1:2n]) # current distribution "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 0.6.0",
   "language": "julia",
   "name": "julia-0.6"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "0.6.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
